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ALTERNATING  DIMENSION  PLASMA  TRANSPORT  IN  THREE  DIMENSIONS 

Harold  Grad 
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Mew  York  University 

New  York,  NY   10012/USA 

ABSTRACT 

The  alternating   dimension    (1^  D)  method  of  solving  macroscopic 
adicbatic  and  transport  problems  is  here  generalized  to  arbitrary  3-D 
'.-.oroidal  plasma  confinement  systems  =   The  principal  new  result  is  the 
derivation  of  an  evolution  equation  for  the  poloidal  and  toroidal  fluxes 
in  which  second  derivatives  can  be  explicitly  exhibited  to  show  that 
the  system  is  diffusive.   This  extends  previous  results  in  2-D,  axial 
symmetry  and  helical  symmetry,  where  the  flux  functions  for  the  magnetic 
field  are  explicit  consequences  of  an  ignorable  coordinate  fand  the  EBT 
closed  magnetic  line  configuration.   The  eigenvalues  {diffusion  coeffi- 
cients) are  evaluated  and  are  shown  to  represent  one-dimensional  relative 
diffusion   among  the  adiabatic  variables,  independent  of  the  representa- 
tion (e.g.  whether  diffusion  is  measured  relative  to  mass,  or  toroidal 
flux,  or  poloidal  flux).   The  skin  effect  diffusion  coefficient  decouples 
from  the  other  coefficients  and  represents  diffusion  of  one  magnetic 
field  component  relative  to  the  other.   Other  transport  coefficients 
such  as  those  for  mass  and  energy  flow  are  intrinsically  coupled.   As 
in  previously  implemented  alternating  dimension  codes,  a  3-D  code  built 
to  these  specifications  should  be  expected  to  be  extremely  accurate  and 
efficient. 

1.   INTRODUCTION 

The  macroscopic  evolution  of  a  plasma  is  governed  by  some  form  of 
conservation  of  mass,  momentum,  energy  and  magnetic  flux.   In  both  ideal 
and  dissipative  dynamic  models,  each  conservation  law  is  an  evolution 
eqi'stion  with  a  time  derivative  of  density,  velocity,  temperature,  or 
magnetic  field.   A  model  for  slow   evolution  [1,2,3]  is  obtained  by 
dropping  inertia,  pdu/dc  =  p(3u/3t  +  u-7u),   from  the  momentum  equation, 


(l.L)  Vp  =  J  x  B,    J  =  curl  B 


Formally,  this  can  be  obtained  by  placing  a  small  parameter  in  front 
of  all  time  derivatives,  velocity  components,  and  transport  coeffici- 
ents; all  equations  are  homogeneous  in  9/9t,  u,  and  transport  coeffi- 
cients and  are  unchanged,  except  for  the  equation  of  momentum  conserva- 
tion which  reduces  to  (1.1) „   The  interpretation  of  the  slow  evolution 
model  is  that  one  is  following  a  time  sequence  of  equilibria,  (1.1). 
Without  transport  (the  adiabatic  model) ,  the  equilibrium  varies  as  a 
consequence  of  changing  external  constraints  or  boundary  conditions. 
Transport  causes  an  additional  slow  evolution  through  static  equili- 
brium states.   The  crucial  point  is  that  the  equilibrium  equation  (1.1) 
is  incomplete  (in  a  toroidal,  closed  magnetic  surface  configuration, 
which  is  all  we  shall  consider)   when  only  boundary  conditions  are 
specified.   In  two  dimensions,  for  example   (n  is  the  unit  vector  in 
the  ignorable  z-direction) , 

(1.2)  B  =  n  x  Vi|i  +  n  B   , 
and  equation  (1.1)  reduces  to 

(1.3)  Aty  =  -  dp/dij;  -  f  df/d* 

where  specifying  the  two  profiles  p  =  p(^)  and  B   =  f  (ijj)  reduces  the 
equilibrium  problem  to  a  nonlinear  elliptic  equation.   The  two  profiles 
p(ip,t),  f(i^,t)  will  vary  with  time  as  a  consequence  of  external  con- 
straints or  dissipation. 

In  the  adiabatic    problem,  the  evolution  of  p(i|>,t)  and  f(\p,t)    will 
be  governed  by  specifying  fixed  adiabatic  profiles,  e.g„   UW  and  v(ij>) 

f   y  =  p/(*')Y 

(1.4)  { 

{       v   =  l/V 

Here   i^'(V)  is  the  derivative  with  respect  to  the  volume  V  (area  in  two 
dimensions),  where  if»(V)  is  the  inverse  function  to  V(iJj),  the  volume 
within  a  surface  \p    =  const.   Eliminating  p  and  f  in  favor  of  \i    and  v 
gives  the  generalized  (or  queer)    differential  equation  (QDE) 

(1.5)  Ai|)  =  F(V,^,4>'  A") 

where  the  right  side  is  a  second  order  ordinary  differential   operator 
in  tp  (V)  .   A  certain  amount  of  analytic  theory  has  been  done  [4,5,6], 
and  there  are  very  effective  numerical  methods  [4,7-14]  for  solving 
the  generalized  differential  equation  (1.5).   For  the  moment,  assume 


that  the  problem  of  finding  an  equilibrium  solution  with  given  adiabatic 
profiles  u(<l>)  and  v(t^)  is  fully  understood.   V7e  now  show  how  to  reduce 
the  problem  of  slow  transport  to  the  adiabatic  problem  (up  to  now  done 
only  in  configurations  with  an  ignorable  coordinate  [15]). 

The  key  to  solving  the  transport  problem  is  to  eliminate  the  veloc- 
ity (appearing  in  the  mass,  energy  and  flux  equations)  which  no  longer 
evolves  in  accordance  with  a  given  value  of  3u/3t  but  is  determined  by 
the  constraint  that  p,  B,  and  J  (which  do  satisfy  evolution  equations) 
must  be  in  pressure  balance  at  all  times,  Vp  =  J  x  B.  The  velocity  is 
eliminated  by  first  taking  in  each  transport  equation  a  microcanonicai 
volume  average  on  a  flux  surface,  and  then  introducing  adiabatic  inde- 
pendent and  dependent  variables.  The  result  is  a  set  of  one-dimensional 
equations  for  the  mutual  diffusion  of  all  adiabatic  variables.  We  may 
take,  for  example,  as  a  complete  set  of  adiabatic  variables  the  extensive 


variables,  mass,   M  = 


adM,   and  two  fluxes,   ij> 


pdv,   entropy  S  = 
(poloidal) ,  and  x  (toroidal)  [note  that  u'and  v  in  (1.4)  are  given  by 
u  =  a  (dM/dtJO Y,  v  =  dx/dtjj   in  terms  of  the  present  adiabatic  set].  Taking 
M  as  the  independent  variable,  S (M) ,  ty (M) ,  and  x (M)  are  invariant  pro- 
files in  an  adiabatic  problem,  and  the  mass  densities  3S/3M,  3^/311, 

3x/3M   satisfy  a  set  of  three  coupled  second  order  diffusion  equations 

('1) 
in  M  and  t.   The  three  eigenvalues,  A.  ',  of  the  matrix  of  coefficients 

2     2       1 
of  second  derivatives  (3  c/3M  ,  etc.)  are  the  diffusion   coeffvcvents , 

2 
with  dimensionality  M  /t.   The  values  of  these  eigenvalues  are  essenti- 
ally independent  of  the  choice  of  independent  variable  (see  Appendix) . 

For  example,  taking  -ty   as  the  independent  variable,   3M/3^,  v,  and  a 

2 
satisfy  a  diffusion  system  with  the  eigenvalues  (dimensionally  ip  /t) 

\W  =  (3*/3M)2a{M)  . 

To  be  more  precise,  the  1-D  diffusion  system  with  is  obtained  by 
averaging  is  not  self-contained  since  it  involves  geometrical  coeffi- 
cients such  as  the  inductance  of  a  shell,  and  thermal  or  electrical 
resistivity  of  a  shell.   Also,  the  three  diffusive  equations  contain  a 
fourth  dependent  variable,   viz.,  the  Jacobian  from  volume  to  the  inde- 
pendent  variable,  e.g.   M' (V)  =  p   or   x'(V),  etc   The  equation  count 
(after  losing  an  equation    by  introducing  an  adiabatic  independent 
variable)  is  maintained  by  supplementing  the  1-D  diffusive  system  by 
the  average  pressure  balance. 

The  common  key  to  solving  either  the  adiabatic  or  the  dissipative 
problem  is  to  alternate  between  solution  of  the  1-D  system  and  deter- 
mination of  the  geometrical  coefficients  in  the  1-D  system  from  the  2-D 


(or  3-D)  equilibrium  (this  is  the  Alternating   Dimension   algorithm ,  at 
one  time  called  "l^-D"  algorithm  [3,4,13]),   For  the  adiabatic  problem 
the  entire  1-D  system  is  the  average  pressure  balance  which  is  a  second 
order  ordinary  differential  equation  for  ty  (V)  given  u(40  and  v(\p)    and 
the  appropriate  inductance  coefficients  L. .(V).   Given  an  approximation 
to  L..(V)  [y  (iM  and  v(\\j)    are  fixed]  the  ordinary  differential  equation 
gives  an  approximation  to  ip  (V)  and  to  the  current  [F(iJj)  in  (1„5)];  this 
is  a   inhomogeneous  term  for  an  elliptic  equation  for  ij;(x,y).   From  such 
an  approximation  to  the  current  distribution  in  physical  space  (2-D  or 
3-D) ,  one  finds  the  flux  surfaces  and  corresponding  inductance  coeffi- 
cients and  repeats  the  procedure. 

This  iterative  procedure  has  been  found  to  converge  much  more 
rapidly  than  standard  elliptic  iteration,  e.g.  than  AiJj   ,  =  F(iJ>  )  for 
given  F(<J;).   In  the  dissipative  problem,  the  2-D  (or  3-D)  geometry  has 
to  be  recomputed  only  after  hundreds  of  diffusion  time  steps.  Thus  the 
diffusion  problem  is  almost  entirely  1-D  with  occasional  corrections 
to  the  geometric  coefficients  from  2-D  or  3-D. 

What  is  done  in  this  paper  is: 

(1)  find  an  evolution  equation  for  the  flux  (3ip/3t,  3x/3t)  in  a 
general,  3-D  configuration; 

(2)  identify  this  evolution  equation  as  diffusive  by  explicitly 
identifying  second  derivatives  of  the  adiabatic  variables; 

(3)  calculate  the  eigenvalues  (diffusion  coefficients)  in  complete 
generality; 

(4)  Show  that  the  diffusion  coefficients  are  independent  of  the  repre- 
sentation; in  other  words,  the  physical  picture  is  1-D  diffusion  of 
adiabatic  variables   relative  to  one  another  (note  that  volume  is  not 

a  proper  variable,   nor  is  mass  flow  relative  to  a  fixed  frame  or  rela- 
tive to  V)  . 

At  least  one  3-D  code  of  this  type  is  being  constructed  [16] »  One 
important  advantage  of  the  Alternating  Dimension  algorithm,  not 
considered  in  this  paper,  is  the  efficient  treatment  of  complex  and 
changing  topology  [17,4,8] , 


Toroidal  Configuration 


2.   PRELIMINARIES 


axis 
A 


Figure  1 


Consider  a  family  of  nested  toroidal  surfaces  with  a  simple  closed 
curve,  A,  as  axis.   Take  C1  and  C0  ,  as  indicated  in  Fig.  1,  as  a  basis 
for"  the  homology  of  one-dimensional  closed  curves  on  a  toroidal  surface, 
oriented  such  that  (C,,C0,n)  is  a  right-handed  system  (n  is  the  outward 
normal  to  the  surface) .   The  open  surfaces  E,  and  E_  are  taken  as  in 
the  figure,  with  oriented  boundaries 


(2.1) 


*z±   -  -  c2  , 


3E2  =  C1  -  A 


The  intersection   numbers  are  given  by   {C.,E.}  =  6... 

Surface  Potentials  and  Feriods  [18] 

Consider  the  family  of  nested  toroidal  surfaces  to   be  parameter- 
ized by  the  values  of  a  parameter  p  through  the  level  surfaces  of  a 
given  function  p(x,y,z)  =  const.   Suppose  X'  is  a  single-valued  two- 
dimensional  vector  field  on  a  surface  p.   We  use  the  notation  X'  =  V'a 
for  the  surface   gradient    of  a  surface   potential    a   defined  by  the  property 


(2.2) 


<)  X'  «dx  =  0 


for  all  closed  curves  C  homologous  to  zero  on  the  surface.  The  surface 
periods 


(2.3) 


[a] .    =   V'a«dx 


Ci 


are  path  independent  on  the  given  surface,  but  are  functions  of  p.   This 
is  in  contrast  to  a  conventional  (three-dimensional)  gradient,  X  =  7<J>. 


Its  projection,  7'<j>,  on  a  surface  is  a  special  case  of  a  surface  grad- 
ient with  periods  which  are  independent  of  p.   Note  that  7<{>  is  single- 
valued  (and   7<j>«dx  path-independent)   within  the  three  dimensional 
domain,  whereas  Va  is  single-valued  only  in  a  domain  cut  by  E,  and  £_, 
or  when  projected  on  each  surface.   For  a  three-dimensional  vector 
field  X,  the  property  n«curl  X  =  0  on  a  p-surface  implies  that  the 
projection  X*  is  a  surface  gradient. 

A  very  useful  identity  involving  two  arbitrary  surface  potentials 
a  and  B  is  (Hodge's  Theorem) 

(2.4)  6   7a  x  VB-dS  =  [a]  ^Q]  2   -  [a]2[B]1 

P 
where  the  integral  is  taken  over  a  toroidal  p-surface. 

Next,  suppose  there  is  given  a  vector  field  B  satisfying 

div  B  =  0 
(2.5) 

B«Vp  =  0 

We  define  the  magnetic  flux  periods 


i-J 


(2.6)  ^.  =    B«dS 

E. 

l 

Clearly,   i|» .  =  <J>.(p)  and,  using  p  =  p(x,y,z),  we  can  extend  <p.  as  point 
functions  ijj.(x,y,z).   The  derivative   diK/dtK  (or  its  reciprocal)  is 
the  rotation   number   of  a  B-line  on  a  given  toroidal  surface  (depending 
on  whether  the  rotation  number  is  measured  by  crossings  of  E,  or  E2) • 
For  each  of  the  tjj.(x,y,z),  surface  potentials  a.  can  be  defined 
(on  each  surface  p)  such  that 

(2.7)  B  =  7a.  x  Vip .  . 

From  the  identity  B«dS  =  da  d'jj,  the  following  period  relations  are 
easily  established: 

f   [a,  ],  =  -diJ;9/dtK  ,   [a,K  =  1, 

(2.8)  {     L    X  Z        L  L 

[       [a2]1  =  -1,  [a2]2  =  d*1/di(;2 

Furthermore,  a,  =  (diK/dij;,)  a2  +  aQ(p)  pointwise  (a1  and  a2  are  unique 
within  an  added  constant  on  each  surface) .   We  shall  sometimes  use  the 
notation  (ij>,x)  f°r  (4s /'J^  or  ^o/^i)  in  which  case  we  shall  write 

(2.9)  B  =  7a  x  7^  =  7B  x  Vx 


(2.10) 


To  continue,  assume  now  that  B  and  p  satisfy 

(Jxd  =  Vp,    J  =  curl  B 
div  B  =  0 


[this  is  compatible  with  (2.5)].   Since  n«curl  B  =  0,  one  can  introduce 
the  surface  potential  <\> , 


(2.11) 


B  =  V 


Defining  the  currents 


(2.12) 


Ii  =    J-d 


and  mmf ' s  (magneto-motive  forces) 


I  . 

l 


(2.13) 


<J>   =  0  B«dx  , 
A    J 

A 


(A  is  the  axis)  we  verify 


=    0    B«dx    = 

C. 

l 


wL 


(2.14) 


I      =   -    <J> 
±1  *2 


I      =    $      -    $ 
A2         *1         *A 


From  div   J   =    0   and   J-Vp   =    0   we   can    introduce 


(2.15) 


J    =    V£    x    Vp 


Using  the  identity  J'«dS  =  d£  dp   and  the  notation  £ .  =  [£].,  we  easily 
verify 


(2.16) 


?1  =  -dl2/dp,       c2  =  dl1/dp 


(2.17) 


T,  .     =  -  d$  .  /dp 
1  1 '      * 


Other  useful  identities  are 


(2.18) 
(2.19) 


J«Va.  =  -dp/dij>. 


B«V£  = 


(2.20) 


3  =  C  Vp  +  7- 


and  the  Jacobian  identities 


(2.21) 


2  -  2 

Vai    x    V4K     •    7$   =   B    ,       dV   =   dot    dty   d$/B 


(2.22) 


Vai    x    Vt^    •    V?    =   1    ,      dV   =   da   d\p   d£ 


To  interpret  these  as  Jacobians  in  three  dimensions,  it  is  simplest  to 
cut  the  interior  of  the  torus  at  I.  and  E_  ,  which  leaves  ($,a,tH  or 
(£,a,t{j)  as  single-valued  coordinates  in  the  cut  domain. 

As  an  example,  using  (a,i^,<J>)  as  variables  we  can  evaluate 


(2.23) 

and  write 


?*  =  VB   , 


J'B  d\p 
^cx  '   B2  dp   ' 


(BxVcfVO/B' 


(2.24) 

Surface  Averages  [18,2,4] 


j  =  d£        B*^ 
dijj   ^"a      a2 


The  torus  '|j  <  const,  has  a  volume  V(<W  which  can  be  inverted  as 
U<  (V)  .  From  tMx,y,z)  we  can  extend  V(^)  to  a  point  function  V(x,y,z). 
We  recall  the  well  known  identities 


(2.25) 


{   °s 


Vij7|    \\i  '  (V)   ' 


I      dS 

(f)Twr  =  l 


leading  to  the  microcanonical  (volume  weighted)  surface  averages, 


<<\>>    = 


dS 


vv 


(2.26) 


=   (!>•$({, 


dS 


d_ 
dV 


<$>    dV 


where  the  volume  integral  is  extended  over  the  interior  of  the  torus 
V  <  const. 

We  introduce  the  standard  notation  for  derivatives  of  general 
point  functions  and  functions  which  are  constant  on  V  =  const., 


(2.27)    |^  4>(x,y#zft),    *t  H  |^  ij,  (v,  t)  , 


P'  =  |^  P(V,t)  . 


Some  easily  derived  formulas  are 


(2.28) 


<div  X>  =  <X«VV> 


=  (div    (*-^)\ 

ivvr  / 


ivv| 

(2.30)  <|^>    =    0 

(2.31)  <|i>    =    <$>t   +    <^   |V> 

<2-^>  Jl  -  f- cSI  -  ♦J 

where  <{>(x,y,z,t)  and  X(x,y,z,t)  represent  arbitrary  scalar  and  vector 

functions  respectively    In  (2.32),  ip  can  be  replaced  by  any  function 

which  is  constant  on  V  surfaces  (e„g0  p)  =  For  f  =  f(v,t),  we  have 
the  important  special  case  of  (2.31), 

(2.33)  <§£>  =  ft 

This  together  v/ith  (2.28)  is  the  reason  why  volume  is  the  natural 
independent  variable  in  averaged  conservation  equations. 

The  formula  (2.4)  for  two  surface  potentials  can  now  be  written 

(2.34)  <Va  *  V3  •  7V>  =  [a]]_[6]2  -  [a]2[3]]_  . 
Applying  this  identity  to  (2„21),  we  obtain 

(2.35)  <B2>  =  ^j  +  $2^2 
Similarly,  from  (2018), 


(2.36)  -  p'  =  $1^1  +  $2^2 
and  from  (2.15)  , 

(2.37)  <J«B>  =  <t20J  -  $^2    =  ^ij  +  $2I2 

For  later  purposes,  it  is  important  to  note  that  there  is  no  simple 

2 

formula  involving  periods  for  <J  >  „ 

Surface  Harmonics  [18,19] 

If  two  surface  potentials  f  and  g  are  related  by 

(2.38)  v'f  =  n  x  v'g 
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where  n  is  the  normal  to  the  surface,  then  f  and  g  are  conjugate 
harmonics    in  the  surface  metric.  Similarly,  defining  oj  =  a,ij>,  ,  from 

(2.39)  B  =  Voj  *  VV  =  V* 

we  recognize  that  B  is  a  generalized  harmonic  vector  in  the  surface 
V  =  const,  with  a)  and  $  as  conjugate  surface  harmonics;  w  and  <J  satisfy 
an  elliptic  equation  which  differs  from  the  Cauchy-Riemann  equations 
(2.38)  by  the  appearance  of  a  weight  function  |VV|  presumed  given.  If 
the  surface  and  the  weight  | VV |  are  given,  then  B  is  uniquely  deter- 
mined by  a  pair  of  periods,  either  [w].  (^'  and  \pl)    or  $..  .  From  this 
follows  the  existence  of  2x2  symmetric,  positive  definite  surface 

inductance  and  susceptance  (or  capacitance)  matrices  L  =  L .   and 

-1  1J 

A  =  L   =  A. .  respectively,  relating  the  conjugate  periods 

i>'.   =  l.  .*  . 

(2.40)  X  1]  3 

$.  =  A.  .\1>\ 

From  (2.35),  the  magnetic  energy  density  per  volume  is 
4  B2>  =  \    (^  +  »2*') 

(2.41)  =7*iLij*j 

The  inductance  coefficients  are  purely  geometrical,  depending 
only  on  the  function  V(x,y,z)  [ | VV j  on  a  given  surface]. 

The  average   pressure  balance  (2.36)  is  most  useful  when  written 
in  terms  of  a  mixed  basis  (iK' ,$_)  or  (<|»2,*-|)  (eliminating  the  other  two 
periods  by  the  inductance  formulas  (2.40)], 

d<£ 

(2-42)    -  %  =  (^An)'  +  *2  aijf  /A22  "  $2  a*j-  (Li2/Lii'  • 

and  the  same  formula  interchanging  indices  1  and  2„   In  configurations 
with  an  ignorable  coordinate,  (2C42)  is  the  average  of  an  elliptic 
equation  for  ijj.  .   For  example  in  two  dimensions  [cf.  (1.3)],  (ibj/L-,,)' 
is  the  average  of  AiJ>,  ,  $_  =  f(»K)  and  A22  =  1  ,   A,2  =  0  . 
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3.   CONSERVATION  OF  MASS  AND  ENERGY 

The  conservation  of  mass,  including  sources,  takes  the  form 

(3.1)  |£  +  div  (pu)  =  Z 

We  assume  that  p  is  constant  on  pressure  surfaces  (this  will  follow 
from  the  fact  that  it  holds  for  pressure   and  temperature  [14]). 
Defining 


(3.2)  U  =  m  u«dS  =  <u«VV> 
We  compute  the  averaged  mass  equation 

(3.3)  pt  +  (pU) '  =   <£> 

Note  that  the  independent  variable  V  is  defined  by  p  =  const   which 
has  not  yet  appeared  in  the  present  discussion.   Defining  total  mass  by 


(3.4)  M  =   pdV,    M1  =  p 


we  obtain 


(3.5)  Mfc  +  Il'U  = 


ZdV 


Without  sources  mass  is  conserved,  and  the  derivative  holding  mass 
fixed,  d/dt,  is 

(3.6)  ^L  E  |-  <MM,t)  =  <j>t  +  u<p' 

We  have 
,,  „  dM    . 

(3-7)  ar  =  ° 

(3.8)  ™   =  U 

Mote  that,  even  without  sources,  (3.1)  dees  not  imply  that  the 
flow  velocity,  u,  carries  p  =  const,  surfaces  into  p  =  const.;  (nor 
does  u  carry  ty   or  p) .   M  is  defined  as  the  mass  within  a  V  surface. 
This  is  not  a  Lagrangian  variable  in  (x,y,z)  since  the  fluid  element 
within  a  V  surface  moving  as  in  (3.8)  [or  in  a  fixed  domain  M(x,y,z,t) 
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<  const.]  is  not  preserved.   However   (with  X  =  0)  M   is  a  Lagrangian 
variable  in  (V,t)„ 

The  pointwise  conservation  of  energy  takes  the  form 

(3.9)  |£  +  u'Vp  +  yp  div  u  =  (y-1)  div(K  VT)  +  (y-DnJ2  +  (y-l)Q 

a  2  2       2   1/2  2 

where  T  =  p/p  is  the  temperature,  k  =  k(T)/cj  t   =  k,  p  /T   B   is  the 

-3/2 
transverse  heat  flow  (k   is  an  absolute  constant),  n  =  r,,T     is  the 

resistivity,  and  Q  is  an  energy  source  (per  volume) .   Averaging,  under 

the  assumption  that  T  =  const,  coincides  with  p  =  const,,  gives 

(3.10)  pt  +  up'  +  ypu'  =  (y-1)(<*t')'  +  (y-1)  <nJ2>  +  (y-1)<Q> 

where 

mil  _kiP2  |vv|2       <iP2  ^\n\\ 

(3.11)  k*  =  -T7I  < r->   =  -Tn <1_JL1_> 

T17^     B^        TX/Z(^')      B^ 

Introducing  the  entropy  per  unit  mass  a  =  Cy   log  (p/pY) ,  yields 

(3-12)  c  (y-D  dT  =  (K*T,) '  +  <nj2>   +  <Q> 

Without  dissipation  or  sources  we  have 

***   =  0     ^-  =  0 
dt   u'    dt   u 

M  and  a   are  adiabatic  variables  and  a (M)  or  M(o)  can  be  taken  as  a  fixed 
adiabatic  profile. 

If  the  motion  is  incompressible,  div  u  =  0  [or  incompressible  on 
the  average,  U(p)  =  0],  then  V  is  also  an  adiabatic  variable „   This  is 
the  reason  why  in  some  models  (e„g.  a  low  3,  small  pcloidal  field, 
straight  Tokamak   [2,8]),  V  has  been  used  successfully  as  the  indepen- 
dent 1-D  variable. 
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4.   FLUX  CONSERVATION 


To  obtain  appropriately  averaged  flux  conservation  equations, 
we  use  Maxwell's  equations 

(4.1)  ~   +  curl  E  =  0  ,    div  B  =  0 

0  t 

and  Ohm's  law 


(4.2)  E  +  u  x  B  =  n  J 

Letting  B  =  Va  *  ViJ;, 


9a  n     d\p 


curl  (||  7*  -  |f  Va) 


and  from  (4.2)  , 


(4.3)  E  =  Va  ||  -  V*  |£  +  Vf  . 

To  begin,  f  must  be  restricted  to  the  cut  domain  (Z.  and  E_) ;  project- 
ing (4.3)  on  a  ty    surface  shows  that  V'f  is  single-valued,  thus  f  is  a 
surface  potential..  Adding 

u'x  B  =  Va(u  •  Vip)  -  V<|/(u  •  Va] 

to  (4.3)  yields 

(4.4)  n  J  =  Va(||-  +  u«Vi^)  -  7*  (|£  +  u«Va)  +  Vf  . 

Equations  (4.3)  and  (4.4)  are  valid  with  (B,x)  replacing  (a,ty)    and  a 
related  potential  g  instead  of  f. 
We  remark  that 


(4.5)  <f!-  +  u-V'4>>  =  \\).     +   Ui^' 


dt 


With  this  as  motivation,  take  the  scalar  product  of  (4.4)  first  with  J, 
then  with  B,  and  then  average,  to  obtain 

<nJ2>  =  -^|+  <J-vf> 
(4.6)  d*  dt 


=  -|l+P,'f1^2V 
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<nJ*B>  =  <B'Vf> 
(4.7) 

=  ip'  (£±0-2   ~    f2al^ 

It  remains  to  evaluate  the  periods  f,  and  f_.   Taking  the  jump  of  (4.3) 
across  either  E,  or  £», 

0  =  [Va]  ||  -  V*  [|£]  +  [Vf]  . 

Note  that  only  the  projection  of  [Va]  is  zero  for  a  surface  potential 
a.  First  take  a  =  a.  and  the  indicated  jump  across  Z~   where  [a.]_=  1 
[cf.  (2.8)].   This  implies  that  [Va,]-  =  0  and  [3a  /at],  =  0,  hence 
[f]y   =   c, (t) .  Similarly,  for  a  =  a~  ,  [g] ,  =  -c-(t).  The  remaining  period 
of  f  or  g  can  then  be  eliminated  between  (4.6)  and  (4.7)  giving 

dip.      dip.     2     dl.  ip!     2     I! 

(4-8)  dt1  =  '  diT  <nJ  >   +  d^  <nJ*B>  "  ci=  "  p^  <nJ  >   +  p^  <nJ-B>-Ci 

These  are  the  two  averaged  flux  conservation  equations.  They 
reduce  to  the  conventional  equations  in  the  case  of  an  ignorable 
coordinate  (two  dimensions,  axial  symmetry,  helical  symmetry),   These 
equations  could  be  used,  in  principle,  to  advance  \p  .    in  time,  since 
the  right  side  is  known  in  terms  of  a  given,  instantaneous  equilibrium 
state,  Vp  =  J*b.   There  is,  however,  no  resemblance  to  a  diffusion 
equation,  since  second  derivatives,  ip'.',  are  not  visible,.   From  (2„14), 
(2.40),  (2.37)  we  see  that  I!,  p1,  and  <J«B>  are  each  linear  expres- 
sios  in  \p"    and  ip".   The  second  term  in   (408)  is  explicit  in  second 

derivatives,  but  it  is  not  at  all  a  conventional  diffusion  operator. 

2 

Furthermore,  <J  >  is  not  expressible  in  terms  of  periods.  Examination 

of  the  special  case  of  axial  symmetry  shows  that  there  is  cancellation 

2 
between  undesirable  terns  arising  from  <J  >  and  <J«B>,  also  that 

geometrical  expressions  ether  than  inductance  arise. 

Before  treating  the  general  case,  note  the  cancellation  in  the 

following  linear  combination  of  the  tv/o  equations  (4.8) 

1   d^2    1   d*l        1  Cl    C2 

(4-9'  jj  at- -  ?j  ar1  -  -  jpr- '"J-B>  +  *i  "  *I 

Also,  from  (2.37)  and  (2„4C) 

(4.10)  <J«B>    =    A(^2>£    -    <i>U>%)     -    A 

where 
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A   =    (*')     (AnA;2    -    A^A^)    +    *>2(A11A22    "    A22A11) 


(4.11) 


and 


+    (^^)2(A12A22    -    A22A[2) 


(4.12)  A    =   A22A11    -    K2±2 

Introducing   iK    as    independent   variable   with   the   notation   D/Dt 
holding   <p.    fixed, 

(4.13)  3<|>/3t    =    D4>/Dt    +    [dtp/dip^)  [dty^dt) 

We   obtain,    for    (4.9) 

2 

D\p2  3    ty2  <p' 

(4'14)  Dt~  =    X0   —2-  +   ^"A   +   C1^-C 

(4.15)  AQ    =   n(^{)2A 

Of  course  this  equation  is  incomplete  without  the  rest  of  the  system, 
but  it  strongly  suggests  that  \\> ~    diffuses  relative  to  ij/,  at  a  rate 

given  by  the  diffusion  coefficient  >..  [similarly,  :j>,  diffuses  reia- 

2 
tive  to  ty-    at  a  rate  given  by  the  diffusion  coefficient  rt(ip')    A], 

Measured  relative  to  V  (this  is  only  dimensional  analysis,  for  there 

is  no  such  diffusion  equation) ,  the  diffusicn  coefficient  is  simply 

nA.   It  is  interesting  to  note  that  A  becomes  large  if  the  flux 

surface  is  rippled,  and  it  may  be  unbounded  at  a  separatrix. 

2   2 
Next  we  turn  to  the  coefficient  <  |  Vij;  |  /B  >  which  appears  in  the 

i   I  2   2      2 
average  heat  conductivity  (3.11)  [in  axial  symmetry  <  j  ViJ;  |  /B  >  =<r  >]. 

2  2 

It  also  makes  an  appearance  as  part  of  <J  >,  viz.  <JX>    (perpendicular 

to  B) 

BxVp 
i=   B2 

<J2>    =    <|Vp|2/B2>=     (^)2<|V^|2/32> 

(4.16)  ||    <J2>       =      E;<|V^|2/B2> 

2   2 

Assuming  that  <  |  ViJ;  |  /B  >  is  geometrical,  the  expression  (4.16),  which 

is  what  appears  in  (4.8),  is  linear  in  second  derivatives  of  \\>    through 
the  factor  d' . 
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We  must  generalize  the  concept  of  geometrical    coefficient  (such 

2 

as  A.  •  ,  or  <r  >  in  axial  symmetry)   to  quasigeometrical      which 

involves  the  geometry  alone  and  first  derivatives,  ty.,    preferably 
as  a  homogeneous  function  of  order  zero  in  tj; !  (i.e„  a  function  of  the 
rotation  number  ip'/ifj1).   Many  properties  of  a  second  order  parabolic 
system  will  be  preserved  with  such  dependence  of  the  coefficients. 

At  first  glance,  the  energy  equation   (3.12)  seems  to  give  second 

2    2 
derivatives  via  (<*)'  ^  < | V^  | /B  >'.   However,  in  accordance  with  the 

adiabatic  format  described  in  the  Introduction  and  in  the  Appendix, 

the  proper   dependent  variables  associated  with  flux  are  ty !  rather  than 

t-;  thus  quasigeometrical   implies  a  coefficient  which  contains  only 

the  dependent  variable  itself,  and  no  derivatives. 

2  2 

To  show  that   <  |  Vip  |  /B  >  is  quasigeometrical  we  recall  that  B 

satisfies  an  elliptic  equation  in  a  surface,   The  manifold  of  solu- 
tions of  this  elliptic  equation  on  a  torus  is  a  two-dimensional 
linear  vector  space,   As  a  basis  in  this  space ^introduce 

B(i)  _  V0(i)  x  vv  =  V<})(1)  ,      i  =  (1,2) 
(4.17) 

[9U']  .  =  6.  . 

3  ID 

From  the  defining  equations  (4.17),  the  B    depend  only  on  the 
geometry,  viz.,  on  the  function  V(x,y,z).   Evidently, 

(4.18)  B  =  i\)'1B{1)    -    ty'2B{2) 

(however,  note  that  these  are  not    Hamada  coordinates).   From  (4.18) 

2   2 
we  conclude  that  <|Vi|j.  |  /B  >   depends  only  on  ty\      as  a  homogeneous 

function  of  order  zero  [and,  of  course,  on  V(x,y,z)]. 

2 

Next  we  turn  to  <J,,>.   On  a  given  surface  with  a  given  vector 

field  B  =  Va  *  V  ip ,   the  current  J  and  its  potential  Z    are  almost 
determined  by  the  ODE  (2.19),  B-7<;  =  1.   To  be  precise  (Ref.  [4], 
Appendix) ,  if  Cg  is  any  solution  of  this  ODE,  the  general  solution 
is  given  by 

(4.19)  x,    =  C,Q    +   cc 

where  c  is  an  arbitrary  constant  (i.e.  function  of  v) .  Similarly, 

(4.20)  J  =  ^o  +  ca|B 

For  a  specific  value  of  c,  c,    in  (4.19)  and  J  in  (4.20)  are  the 
originally  given  equilibrium  quantities.   Note  that  J±    is  independent 
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of  c,  and  JXB  =  JnxB  =  Vp.   There  are  tv/o  choices  of  c  such  that 
one  of  the  periods  C .  (or  <t>  •  )  is  zero;  such  a  choice  has  the  property 
that  the  corresponding  direction  field,  J,  consists  of  closed  lines.. 
These  two  distinguished  values  of  c  lead  to  Cf-»  and  Jf-%  , 

C(1)  =  C  -  [Cl  ^  , 
[c(1)]  =  (l/*{,0) 


(4.21) 
and 


(4.22) 


^(2)  =  ?  +  ^jl0t2 


d$. 
J(l)  =  J  +  a^ 


I(1)  =  (0,  -pV*[), 


[?(2)1  =  (0,  1/<|»2) 
d*. 


:2> 


=  j  - 


d^2 

1(2)  ■  (PV^2,0) 


(this  representation  is    related  to  Hamada  coordinates) .  From  the 
construction  of  these  solutions  of  the  ODE,  3*7^  =  1,  with  one  period 
of  c  zero,  it  is  clear  that  the   potentials  £,   depend  only  on  the 
geometry  and  the  periods  i|>!  of  B.   Consistent  with  this,  the  periods 
of  ?  / -I*  and  C/-jv  also  depend  only  on  !(/!  (4.21),  whereas  the  periods 
of  c,    itself,  l!/p'  contain   second  derivatives.   In  other  words, 
the  second  derivatives  enter  into  solution  of  the  ODE  only  as  speci- 
fied side  conditions*  Since 


(i)  =  vC(i)  x  Vp  =  p.  v5(i)  x  VV, 


J.../p'  is  quasigeometric  and  homogeneous  of  order  -1  in  ij/!. 

This  is  sufficient  to  evaluate  <J  >  in  terms  of  <J,..>  and 
periods;  after  some  manipulation, 

(4.23) 

(4.24) 


l^i 
n    dt 


dip  -  2  i        ci 

3— -  <J,.,>    -    <J>  .  <f> ! /ij> . (j    ?   i) 

dp  d)  f]   ]/ri     n 


d<t>  . 
3 


=    -    *E_   R      -    . 

dt^i      i  3    d\l>L 


(j    t    i) 


where  R.  is  quasigeometric  and  can  be  written  as 


(4.25) 


R.  =  S.  +  T. 
ill 

S.  =  <  |  Vip-   I  2/B2> 
1      i  f1 1  / 


•   =  <(|5_  _  r  )V>  ,    T,  =  <(|^-  +  C.)V 

1     *-3a,     2-1  2     vda,     1J 


The  second  derivatives  in  (4.24)  are  now  explicitly  visible  in 
dp/dt/;.  and  dp./d;..   For  later  use  we  introduce 


(4.26)  /  X  I         - 
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R,    =  R.    -    <$>'i/<B2>  (j    /   i) 


T.     =   T.    -    (J)V<B    >  (j    /    i) 


and   note    that 


(4.27)  Ri    >    Ti    >    0 

(4.28)  (*2)2TX    =    (^{)2T2 


The    inequality   T.     >    0    follows    from 

?    ?         2  2    2 

<f    B    ><B    >       >       <f    B    > 


and  the  identity  is  merely  tedious. 


5.   COMPLETE  ONE-DIMENSIONAL  SYSTEM  AND  DIFFUSION  COEFFICIENTS 

In  terms  of  V  as  independent  variable,  the  complete  system  of 
averaged  equations,  collected  from  (2.3),  (3.12),  (4.23),  and  (2.36) 
is  (dropping  £  and  Q) 

(5.1)  p      +    (pU)?    =0,         Mfc   +   UM'    =    0 

(5.2)  c    (P_1)     (at    +    Uo')    =    U*T')«    +    n<J2> 

di|>. 

(5.3)  __L=    _    _IL,    (p.R.     +    *.„,.,     -c.  (j    /    i) 

(5.4)  ~  P*  =  $[$[    +    ♦2*2 
Since  p  is  not  an  adiabatic  variable,  introduce 

(5.5)  q  =  o/i>[   =  314/3*! 
and  combine  (5.1)  with  (5.3)  to  obtain 

(5'6)  -q^-i^Bx  <*!?'  +*j*j)}' 

Only   one  of  the  two  equations  (5.6)  is  to  be  used  to  replace  (5.1) 
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if  \p  .    or 


ad'A      is  taken  as  the  independent  variable;  both  are  used 
[replacing  (5.3)]   if  M  is  the  independent  variable. 

To  be  specific,  choose  M  as  the  independent  variable.   Using  the 
notation 

(5.7)  $  a  9<J)/8ri  =  (fj'/P 

(5.8)  i>L=    f±  =  l/5± 

we  obtain  for  the  adiabatic  variables  (a,f.)  and  the  auxiliary  vari- 
able p , 

(5.9)  c^ITdT  =  p(k*p*>*  +  n<J2> 

(5.10)  —^  =  -  (nRiP/fi)  -  (ncf^j/fi)     (j  ?   i) 

(5.11)  -  p  =  ptf^  +  f2*2) 

with  the  supplementary  relations 

p   =    pY   exp    (a/Cv) 

Y-l 

T     =   p'        exp    (a/c    ) 

*.=  p  y  a. .f . 

<B2>   =      p(f14>1   +   £2*2) 

<J«B>  =   p(4>2,ll1  "  ^i1^ 
<J  >   =  (I^/i//  )  p     +  <J.B>2/<B2>  ,   (i=l,2) 

2 
There  is  a  choice  of  two  expressions  for  <J  >.   Noting 

(5.12)  p/p  =  y  P/P  +  ^/cv 

we  see  that  (5.11)  expresses  p  as  a  linear  combination  of  (a,f.) 
which  allows  p  to  be  eliminated  from  (5.9)  and  (5.10)  in  favor  of 
(o,f.)   which  are  now  explicit. 

A  tedious  calculation  yields  the  eigenvalues  of  the  3x3  matrix 
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of  coefficients  of  second  derivatives,  AQ  and  A  +  , 

(5.13)  AQ  =  hp2a 

1/2 

(5.14)  A±  =  i  (Ax  +    X2)    +  i  [xj  +  A2,  -  26A1A2] 

where 

(5.15)  6    =   P+(2/Y  -  ^)<B2>  ,         0<6<1 

p+<B  > 

2 

(5.16)  A,  =  n(R./f2)  YP<B  p 

x       -1   1   YP+<B  > 

2 

(5.17)  A-  =  (Y-l)K,(S./f2)  -4to   P+<V 

*  x   l   i   Tx/^  Yp+<B^> 

A. 

The  geometric  factors  R.  and  S.  are  both  positive   [cf.  (4.27)];  the 

combinations  R./f   /  S./ff  are  independent  of  i   [cf.  (4.28)].  These 

i   i     i   i 
eigenvalues  agree  with  previously  calculated  special  cases  in  axial 

symmetry  [15]  and  for  EBT  [20] .   There  is  a  more  sophisticated  trans- 
port model  (tensor  resistivity  and  heat  conductivity)  which  has  been 
carried  out  in  axial  symmetry  and  which  involves  additional  geometri- 
cal coefficients  [14]  ;   so  far  we  have    not   extended   this 
analysis  to  the  general  3-D  geometry. 

The  expressions  A.  and  A_  are  the  contributions  from  resistivity 

and  heat  flow  respectively,  but  they  are  not  the  eigenvalues  (unless 

2 
6  =  1,  which  is  approached  only  at  high  3,  <B  >/p  ■>  0).   The  geometric 

factors  in  A.,  and  A_  approach  their  "Pf  irsch-Schliiter"  values  in  the 

limit  of  a  large  aspect  circular  cross-section  Tokamak.   We  see  that 

these  factors  have  been  recovered  in  the  eigenvalues,  although  they 

were  lost  in  the  mass  flow  and  energy  flow  (which  are  nonlocal)  in  the 

transition  from  Pf irsch-Schluter  to  Grad-Hogan. 
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6.   DISCUSSION 

The  claim  that  we  have  given  a  general  3-D  version  of  classical 
transport  requires  some  interpretation.   It  was  pointed  out  in  1967 
[21,22]  that  the  equation  Vp  =  JxB  is  ill-posed  in  any  toroidal  config- 
uration other  than  2-D,  axial  symmetry,  helical  symmetry,  and  reflection 
symmetry  (EBT) .   The  (now  pervasive)  appearance  of  island  structures 
was  predicted,  and,  for  example,  the  so-called  major   disruption    in  a 
Tokamak  is  evidently  a  disappearance  of  equilibrium  rather  than  an 
instability.   The  problem  of  islation  and  complex  topologies  has  been 
extensively  studied  using  alternating  dimension  codes,  but  so  far  in 
configurations  where  existence  of  solutions  is  not  in  doubt. 

The  practical  use  of  the  present  formulation  lies  in  the  fact 
that  the  failure  of  equilibrium  to  exist  is  sometimes  unimportant  and 
sometimes  invisible.   The  latter  is  a   property  of  some  analytic  expan- 
sions, for  example,  those  in  which  the  rotation  number  per  period  is 
small  (all  resonances  are  invisible  to  such  expansions)  [23,24]. 
Insofar  as  such  expansions  may  approximate  reality,  the  equilibrium 
calculation  can  be  extended  to  include  transport. 

A  second  application  is  to  numerical  equilibrium  calculations  in 
which  the  existence  of  simple  flux  surfaces  is  forced  [25]  .   Refining 
the  mesh  will,  of  course,  lead  to  difficulty;  but  these  difficulties 
may  be  hard  to  identify  in  a  3-D  computation!   Again,  insofar  as  a 
numerical   code  is  a  black  box  which  seems  to  give  numerical  equilibria, 
it  can  be  extended  to  include  transport. 

From  experience  in  2-D,  we  should  expect  great  accuracy  and 
efficiency  for  these  transport  techniques.   As  one  example,  a  computa- 
tion simulating  the  resistive  transfer  from  a  Belt  Pinch  (simple 
topology)  to  a  fully  developed  Doublet  (figure  eight  magnetic  surfaces) , 
with  greater  accuracy  than  is  physically  meaningful,  was  run  in  approxi- 
mately fifteen  seconds  of  7600  time  [8] .   In  another  example,  the  power 
output, 


J 


2 

r,J    (in  this  case  J  is  a  second  derivative  since  ty    is  the 


dependent  variable) ,  was  calculated  for  an  oscillating  separatrix  with 

4 
an  accuracy  of  one  part  in  10   using  a  32x64  2-D  mesh  [26], 

One  early  operating  version  of  the  algorithm   ([17],  1974)  used 

piecewise  constant   (in  time)  geometrical  coefficients.   Later  versions 

use  linear  (in  time)  interpolations  of  geometrical  coefficients  between 

successive  calculations  of  equilibria,  together  with  iteration  of  the 

diffusion  step  until  the  advance  geometry  converges.   To  optimize  a 

3-D  transport  calculation,  one  should  be  prepared  to  minimize  the  number 

of  geometrical  computations  by  using  more  sophisticated  time  dependence. 

However,  in  the  above-quoted  transfer  from  Belt  Pinch  to  Doublet,  using 
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only  five  geometry  time  steps  was   sufficient  to  give  adequate 
accuracy.   In  any  less  radical  change  of  geometry,  two  or  three 
equilibrium  calculations  (with  the  remainder  of  the  computation  in  1-D) 
should  suffice  for  practical  purposes  (but  much  greater  accuracy  is 
required  for  mathematical  purposes  such  as  the  study  of  singularities 
and  boundary  layers) „ 

APPENDIX:   INVARIANCE  OF  DIFFUSION  COEFFICIENTS 


Consider  a  set  of  dependent  variables  u  =  (u, ,...,u  )  which 
satisfy  the  diffusion  system 


(A.l) 


9u.    _        9u. 

*  =  |_  (A.  .  ~-l)     +  R. 


a2u. 

=  Bij  Tx*  +   Ri 


A  and  R  are  functions  of  (x,u,9u/9x);  more  generally,  R  can  also 


include  "lower  order"  terms  such  as 


f(x,u,9u/9x)  dx.   The  principal 


part  (i.e.  excluding  R)  is  in  conservation  form;  the  diffusion  coeffi- 
cients are  defined  to  be  the  eigenvalues  A. (x,u,9u/9x)  of  the  matrix 
B  and  are  assumed  to  be  negative  (B  may  incorporate  second  derivatives 
from  the  expansion  of  9A/9x,  but  in  the  actual  model  in  the  text,  A 
does  not  depend  on  9u/9x) . 


Theorem. 


If  x  is  replaced  by  the  new  independent  variable  £  =  u .dx, 


then  u(£,t)  satisfies  a  system  similar  to  (A.l);  specifically  the  new 

.2,   ....        .,..    _,  „.       ,  =   2, 

i 
We  use  the  single  expression  "R"  to  represent  any  combination  of 

lower  order  terms.   With  the  notation  9<j>  (£,  t) /9t  =  d<j>/dt  , 


principal  part  is  A1  =  u.A;  the  eigenvalues  of  B'  are  A.   :  u"A . . 


dt 


il  21.  li. 

9t  ~  9£  3t 


by  use  of 


21  +   R 
9t  +  R 
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3t 


9u. 


dx 


oU  . 

3 


lj  9x 


R-,  dx 


Also 


9u 
9x 


=  u 


oU 

1  95 


therefore 
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3  (.    3U-,       3  (        A  3u^    S   ,  2  „  3u-> 
33T  (A  3^  =  ul  H  (ul  A  3T}  =  31  ^1  A  3T} 


+  R 


and 


du    3   .  2  .  3u, 

3t  =  3T  (ui   A  3T>  +  R 


A  second  transformation,  replacing  u,  by  1/u,  ,  leaves  the  eigen- 
values unchanged  (this  result  is  classical  for  any  transformation  of 
dependent  variables) ,   We  use  the  two  transformations  on  a  set  of  exten- 
sive adiabatic  variables  (aQ  ,a..  ,  .  .  u  ,  a  )   which  are  assumed  to  satisfy 
(A.l)  where 

x  =  o   ,    u  =  (3a, /3x,  ...,  3a  /3x) 

The  twice-transformed  system  is  equivalent  to  interchanging  an  and  c,  , 


£,      =   a1    =        ux  dx 


u,  = 


1    35      ux  * 

In  other  words,  choosing  any  a.  as  independent  variable  and  3o./3o-, 
j  ^    i,   as  the  set  of  dependent  variables   gives  a  system  with  the 
same  eigenvalues,  independent  of  i,  except  for  a  uniform  factor. 
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